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Image Segmentation and Restoration Using Parametric Contours With 

Free Endpoints 

Heike Benninghoff* and Harald Garcke^ 


Abstract 

In this paper, we introduce a novel approach for active 
contours with free endpoints. A scheme is presented for 
image segmentation and restoration based on a discrete 
version of the Mumford-Shah functional where the con¬ 
tours can be both closed and open curves. Additional to 
a flow of the curves in normal direction, evolution laws for 
the tangential flow of the endpoints are derived. Using 
a parametric approach to describe the evolving contours 
together with an edge-preserving denoising, we obtain a 
fast method for image segmentation and restoration. The 
analytical and numerical schemes are presented followed 
by numerical experiments with artificial test images and 
with a real medical image. 

Keywords: Image segmentation, image restoration, 
active contours, Mumford-Shah, Chan-Vese, paramet¬ 
ric method, variational methods, free endpoints, open 
boundaries. 


1 Introduction 

This article addresses important classical problems in im¬ 
age processing: image segmentation, edge detection and 
image restoration. 

Image segmentation aims at partitioning a given image 
into its constituent parts, also called regions or phases. 
A segmentation of an image can be given by a set of 
region boundaries and edges. Different types of edges 
can occur in images: edges can be boundaries of objects 
and separate these objects from their background or from 
each other. But edges can also end inside the image at a 
location where no other edge continues. 

Boundaries of objects can be modeled with so-called 
interface curves. Non-interface curves are curves which 
do not separate two different regions in the image. Such 
curves have one or two so-called free endpoints. 

Image restoration aims at reducing or removing noise 
which affects a given image. Typically, a blurring of 
the sharp edges in the image should be prevented when 
smoothing an image. This results in the need of an edge 
preserving image denoising method. 
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Image segmentation including edge detection can be 
performed with active contours (also called snakes), first 
proposed by Kass, Witkin, and Terzopoulos |22 in 1988. 
Since this time, the popular method is applied and fur¬ 
ther developed by many authors, e.g. mmmmm 
25plp3 . Using active contours, a curve evolves in order 


to minimize a given energy functional. The energy func¬ 
tional should be designed such that a minimizing curve 
matches with the region boundaries or edges in the image. 

The Mumford-Shah functional 28 can be used for 
both image segmentation and image restoration. A pair 
(T, u) should be found which minimizes the Mumford- 
Shah energy, where T is a set of curves and u is a piece- 
wise smooth function with possible discontinuities across 
T. Having found a solution (T,?/), a segmentation of the 
image is given by the set of object boundaries and edges 
T, and a denoised version of the image is given by the 
piecewise smooth approximation u. 

An important variant of the Mumford-Shah problem 
is the restriction to piecewise constant image approxi¬ 


mations u , the so-called minimal partition problem 15 


However if edges with free endpoints, also called crack- 
tips 28 , occur, the piecewise constant approximation 
will not be applicable. 

It is also possible to approximate the Mumford-Shah 
functional by a sequence of simpler elliptic variational 
problems as introduced by Ambrosio and Tortorelli |l]. 
They replaced the curve T by a 2D function for which a 
phase field type energy is added to the functional. 

Image segmentation and restoration are classical areas 


in image analysis, see 22, 

283334 

, but still significant 

in more present research, 

see e.g. 

3,8 

12 13 15 

17 20 


27||37||38] to mention some selected works. There is also 
a variety of related image processing tasks like object 
detection JH0 or pattern recognition [lO ), feature ex¬ 
traction 


29 and anomaly detection 16 


The image segmentation method, considered and de¬ 
veloped in this article, also uses the evolution of curves. 
The resulting evolution equations, derived from the 
Mumford-Shah functional, can be written as parabolic 
partial differential equations for a parametrization of the 
curves T. The restoration is performed by solving a diffu¬ 
sion equation for u , also derived from the Mumford-Shah 
model. By using the location of the curves T, we obtain 
an edge-preserving smoothing. 

Open active contours, i.e. active contours with free 
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endpoints, are also considered by 


24 


where the au¬ 


thors propose a method for detection of open boundaries 
based on an edge detector which uses the image gradient. 
Here, we consider approaches based on the Mumford- 
Shah model. Using convex relaxation approaches, global 
minimizers of the Mumford-Shah functional are deter¬ 
mined in 32 . The method can also handle free end¬ 


points. In 36 , the level set method is used for evolving 
curves with free endpoints. However, two level set func¬ 
tions and artificial regions are needed to describe a curve 
with free endpoints. 

During the evolution of curves, topology changes like 
splitting or merging can occur, since the number and the 
topology type of edges and region boundaries is often not 
known in advance. Using indirect methods like level set 
and phase field techniques, topology changes are handled 
automatically. It is often argued that the inability to 
change the topology of curves is the main disadvantage 


of parametric methods like the original snake model 22 


In this paper, we extend an efficient method to detect 
and perform topology changes (presented in [9] and based 
on the original idea of [5,26]), such that also topology 
changes of curves with free endpoints can be handled. 

The objective of this article is to solve the Mumford- 
Shah problem including curves with free endpoints with 
a parametric approach. The method we propose is based 
on a parametrization of the evolving curves. We show 
how a method developed for interface curves [9] can be 
extended for curves with free endpoints. With the pre¬ 
sented concept for image segmentation and restoration, 
we can easily process images with both open and closed 
edges. Our method is very efficient from a computational 
point of view, since the curve evolution problem is a one¬ 
dimensional problem and no artificial regions have to be 


used compared to 36 


2 Image Processing with Para¬ 
metric Contours 


Let uo : D —>• M be an image function describing for each 
point in the image domain Del 2 the intensity of the 
image. 


The Mumford-Shah method 28 for optimal approx¬ 


imation of images aims at finding a set of curves T = 
Ti U ... Ttvc and a piecewise smooth function u : fi M 
with possible discontinuities across T approximating the 
original image uq. The energy to be minimized is 

E(u, r) = a\T\ + [ || Vu|| 2 da; + A [ (u 0 - u ) 2 dr, (1) 
Jn \r Jn 


where cr, A > 0 are weighting parameters and \T\ denotes 
the total length of the curves in T. 

A minimizer of the Mumford-Shah functional provides 
(i) a restoration of the possible noisy original image by a 



Figure 1: Image containing an edge with a free endpoint. 


piecewise smooth approximation u and (ii) a segmenta¬ 
tion of the image given by a union of curves T represent¬ 
ing the set of edges in the image. The curves belonging 
to T can be sharp edges where the image function rapidly 
changes, but they can also be so-called weak edges where 


the image function smoothly changes its value, see 15 . 

The contours T^, i = 1,..., Nc , may be closed contours 
with dT = 0, or open contours with two endpoints. The 
endpoints may he on the image boundary <9fi, may be¬ 
long to triple junctions where three curves meet, or may 
be free endpoints, cf. the conjecture of Mumford and 
Shah 128]. In the latter case, the endpoint is a point in¬ 
side the image domain, where no other curve continues. 
Figure [l] shows an image where an edge occurs which 
terminates near the image center. The edge can be rep¬ 
resented by a curve with one endpoint located at the left 
image boundary dfi and one endpoint being a free end¬ 
point, located close to the image center. 

In [9] , we proposed a parametric method for image seg¬ 
mentation with piecewise constant image approximations 
u and interface curves Ti,..., Tn c , each separating two 
regions. There, we considered a decomposition of the im¬ 
age in Nr regions fii,..., fi^n separated by curves 
i = 1,..., iVo, and approximations u\Q k = c& G M. In 
that case, the functional 0 reduces to 


Nr r 

E(T,ci,.. .,c N r ) s= <r|r| + / (u 0 -c k ) 2 dx 

k= 1 


( 2 ) 


15 . Using methods from the theory of calculus of 


variations, the following evolution equation can be de¬ 
rived for time-dependent curves r^(t), t G [0, T\: 


(V n )i = (TKi + Fi, i = 1,..., N c , 


(3) 


where (V n )i is the normal velocity of TV), /%,; is the cur- 
vature, and Fi is given by 

= A[(u 0 - c k+{i) {t)) 2 - (no - c k - {i) (t)) 2 ], (4) 

with t G [0, T\. The indices k + (i),k~ (i) G {1, ...,7V#} 
denote the two regions which are separated by Ti(t). The 
coefficients c&(£), k = 1 ,... ,Nr, are the mean of uo in 
the region Q&(t). 

In practice, the segmentation problem can be solved 
in a two-step approach. For discrete time steps t > 0, 


2 














the coefficients c& are computed using the current set 
of curves. This is followed by an update of the curves 
T(£) T(t + At), performed by solving the evolution 
equation ©• 

For some images, the piecewise constant approximation 
is not applicable, see the exemplary image in Figure [l] 
For such images, the image domain cannot be decom¬ 
posed in regions separated by interface curves. 

Consequently, we modify the two-step approach of |9 , 
such that also non-interface curves with free endpoints 
can be dealt with. In the first step, we will solve a 
diffusion equation in the image domain resulting in a 
piecewise smooth approximation u. Instead of using the 
coefficients c&, we will consider for p E Ti(t) the limit 
u ± (p) = lim e ^ 0 ,e>o u{p db eVi(p )), where Vi is a normal 
vector field on Ti(t). Having computed u , we solve the 
evolution equation © with a modified external term i^, 
using instead of constants c^±^y 

Before, presenting further details, we first consider the 
regularity of a solution of the Mumford-Shah problem 
at the free endpoint. Fixing the curves T, let u denote 
the minimizer of the Mumford-Shah energy ©. At free 
endpoints, problems concerning the regularity of u occur, 
cf. 28 . Expressed in polar coordinates (r, </>) centered at 
the free endpoint, the solution u is of the form 


u(r, <j>) = cr 1/2 sin(I(0 - <j> 0 )) + v(r, 0), (5) 

where v is a C ,1 -function and c,<j >o are constants, see [I]. 

For image segmentation, we later need to solve the 
problem on a discrete set: Let fl h be a rectangular grid 
of nodes covering Q with grid size h > 0. We replace the 
second integral on the right hand side of 0 by a sum 
containing difference quotients of the form 


V % h u{z) = -(u(z + hei) - u(z)), zeQ h , ( 6 ) 


where E M 2 are the standard basis vectors of M 2 , i = 
1,2. For image segmentation applications, we choose the 
pixel grid, i.e. we use h = 1. For the approximating 
sum, we have to exclude terms where the line [£, z-j- hej\ 
intersects with the curve T. 

Instead of the original Mumford-Shah functional 0 , 
we thus consider the energy 


E h (T,u) = cr|r| + £ (1 -a x (z))(V 2 h u(Z)) 2 + 

zen h 

s.t. z+he2€.Q h 

+ E ( 1 -%( 5 '))(V ^(^)) 2 + A f(u 0 -u) 2 dx, 
ten" Jn 

s.t. 

( 7 ) 


2.1 Example 

We consider one single open curve T. Let x : [0,1] —)> M 2 
with £([ 0 , 1 ]) = T be a parameterization of the curve. 
Let x(0) be a free endpoint and let x(l) intersect with 
the image boundary. 

Figure [ 2 ] visualizes a possible situation near the free 

endpoint £(0). Let z ++ , _, z _, £_ + denote the four 

grid points around £(0) as shown in Figure [ 2 j In this 
example, the tangential vector of the curve at £( 0 ) is 
r( 0 ) = £ s ( 0 ) = ei, where s denotes the arc-length of the 
curve. 

Considering z = £4 _, the line [£4 _, £++] and T inter¬ 
sect. Thus, a x (z-\ _) is set to 1. For z — we define 

a factor 


a x (f__):=l((^)i-(x(0))i), (8) 

where (. )i denotes the i-th component of a vector, i = 
1,2. The factor a x (z--) describes how far the curve has 

entered the square given by £++, z. *4_, z _, £_ + . 

For z E £ 7 ^ £_, we set a x (z) = 0 if [£, £+ he 2 ] D 

T = 0 and a x (z) = 1 else. The factor a y (z) is defined 
similarly. 

We now want to vary T in direction —r(0) at £(0). We 
consider a second curve T e , e > 0 , with a parameteriza¬ 
tion £ e , such that £ e (0) = £(0) — er(0). We can assume 
that e is small enough, such that £ e ( 0 ) is still inside the 

square given by £ ++ , £^_, £_, £_ + . 

The energy difference is 

E h (r e , u ) - E h { r, u) =ae + (1 - «,(£__) - e)(V^(5*__)) 2 
- (1 - a x (z__))(V 2 h u(z__)) 2 
=ae-e(V 2 u(z .__)) 2 . 


The energy will decrease if 

a < (V^(z__)) 2 . (9) 

Thus a motion of a curve in direction —r(0) = — e\ at 
the free endpoint £( 0 ) requires that the square of the 
difference quotient of u in ^-direction at £_is suffi¬ 

cient large compared to the weighting parameter a of the 
length term in the energy 0 . 

2.2 General Case 

We consider a curve with one or two free endpoints 
with tangential vector r(p) at the free endpoint £(p), 
p E {0,1}. We define the factors a x and a y as follows: 

Let £ ++ (p),£^_ (p),£ _(p),£_ + (p) denote the four grid 

points around x{p). 

If p = 0 and t(0) . e\ > 0 or if p = 1 and r( 1 ). e i < 0, 
we define z p V := £__ (p) and set 


where a x (z), a y {z ) E [0,1] are scalar terms. If [£, z + he i] 
intersects with T, a y (z) is set to 1 . 
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Figure 2: Illustration of the pixel grid close to the free 
endpoint. 


If p = 0 and t( 0) . e\ < 0 or if p = 1 and r( 1). e\ > 0, 
we define £ P,:L := £+_(p) and set 

MF’ 1 ) := 1 - \ ((F 1 )! - (f(p))0 . 

If p = 0 and r(0). e *2 > 0 or if p = 1 and r( 1). £2 < 0, 
we define z p ,2 := £_(p) and set 

a,(F 2 ) :=l-t((x(p)) 2 -(F 2 ) 2 ). 

If p = 0 and t( 0) . e *2 < 0 or if p = 1 and r( 1). £2 > 0, 
we define z p ’ 2 := z_ + (p) and set 

<*2/(F 2 ) := 1 - t ((F 2 ) 2 - (f(p)) 2 ) . 


Using these definitions, we define the following factors 
for 2 e 


a x {z) 


1, if [£,2 + he 2 \ nr^0, and 

f ^ F 1 , p e {0,1}, 
o^F’ 1 ), if 2 = F\p e {0,1}, 

0, else. 


and 


a y (z) 


1, if [;?,£ +ftei] fl T ^ 0, and 

F 2 , p e {o,i}, 

aty{z p ’ 2 ), if z = Z p ’ 2 ,p e {0,1}, 

0, else. 


For minimizing ([7]), we propose the following approach: 

Assume the case of one curve T with two free endpoints 
parameterized by x : [0,1] M 2 . In the first step, we fix 

u in 0 and consider for fj : [0,1] —> R 2 a variation of x 
of the form x + efj, e > 0. 

Let T 6 " denote the image of x+eff. We use the notation 


E h (ff) := E h (T e,ri ,u) and compute 

E h (T e,ri ,u) 

e=0 

= — a J x ss .ffds — j Fv.fjds+ 

+ ax s (l). ff( 1) - <ra; s (0). ff(0) 

- sign(r(l). ei)rf(l). e^V^u}#’ 1 )) 2 

- sign(f(l). e 2 )ff( 1). e 2 (V^M(£ 1,2 )) 2 
+ sign(f(0). ei)jy(0). ^(V^z 0,1 )) 2 
+ sign(r(0). e 2 )ff( 0). e 2 (V£u(i°’ 2 )) 2 . 

For this computation, integration by parts and a trans¬ 
port theorem are applied. Further, V is a normal vector 
field on T such that the pair (x s , v) is a positive oriented 
basis of M 2 and F is defined as the jump 

F = X[(u 0 - u+) 2 - (u 0 - u~) 2 ]. (10) 



We define the following inner product for functions 
V, X ■ [0,1] ->■ M 2 : 

(*7, x) 2 , r ,ar ■= -X^s + ffil) .\(1)+ff(0)-x(0). (11) 

Now, we consider a family of curves T(t), t G [0, T\. Let 
x : [0,1] x [0, T] —» M 2 be a mapping such that x {., t) is a 
parameterization of T(t), t G [0, T\. We call x a solution 
of the gradient flow equation, if 


Ou,j?) 2 , 


,r,ar 


d 

de 


E h (fJ) 

e=0 


( 12 ) 


holds for all ff : [0,1] M 2 . 

In the following, we consider particular choices of func¬ 
tions ff and derive evolution equations for the curve. 
First, we consider fj = 770 V for a scalar function p 0 : 
[0,1] M with 770 ( 0 ) = 770 ( 1 ) = 0. This provides 



x t .Vr\^ d s 



ax ss .if-\-F)r]Q ds. 


Since 770 is arbitrary chosen (with value 0 at the end¬ 
points), we conclude the following equation for the nor¬ 
mal velocity of the curve: 


V n := xt • v = ctk + F, (13a) 

using the identity 

ku = x ss . (13b) 

Next, we choose ff = 770T, where 770 : [0, 1 ] —>> M is a 
scalar function with 770(0) 7 ^ 0 and 770(1) = 0, i.e. ff( 0) = 
r (0)770(0) and ff( 1 ) = 0. Inserting ff in ([l2|, and using 
p3a> , pbl ) and r = x s , leads to 

xt(0) ■ f( 0 )r? 0 ( 0 ) = cr%( 0 ) 

- sign(f( 0 ). ei)f( 0 )? 7 o( 0 ). e 1 (X7 2 h u(z°’ 1 )) 2 + 

- sign(r(0). e 2 )r(0)rj o (0). e 2 (Vj l u(^ ) ’ 2 )) 2 . 
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Since 770 (0) is arbitrary and sign(r(0). e^)r(0). e* = 
|r(0). ei\ for 7 = 1, 2, we conclude for the tangential ve¬ 
locity 


Vtan(0) :=f t (0) .r(0) 

=cr- |r(0).ei|(V^u(2 0 ’ 1 )) 2 
»|f(0).e 2 |(V^ 2 )) 2 . (13c) 

Choosing 770 (0) = 0 and 770 ( 1 ) 7 ^ 0, we can derive the 
following equation for the tangential velocity in x(l): 

Vtan(l) :=X t (l) ,f(l) 

= -a+|f(l).ei|(V 2 U (f 1 ’ 1 )) 2 
+ 1^(1 )-e 2 |(Vfcw(i' 1 ’ 2 )) 2 . (13d) 


coupled. The cases with triple junctions and boundary 
intersection points are described in |§] in detail. 

We alternately solve the scheme of evolution equations 
(13) and recompute the approximating function u using 
the updated curve set. The function u is obtained by 
solving a diffusion equation on Q h . We note that 0 is 
formulated for a discrete set Q h . We will describe in the 
next section, how u is computed numerically. 


3 Numerical Approximation 

3.1 Numerical Solution of the Evolution 
Equations 


Similarly, choosing fj = 7/0 7, provides the following 
equations for the normal velocity at the free endpoints: 

V„(0) :=x t (0). *7(0) 

= - sign(r(0). elHO). e 1 (V 2 u(i 0 ’ 1 )) 2 

- sign(r(0). e 2 )7(0). e 2 (V^u(5°’ 2 )) 2 , (13e) 


K(l) :=z t (0) .7(0) 

= + sign(f(l). ei)7(l). eKV 2 ^ 1 ’ 1 )) 2 
+ sign(r(l). e 2 )7(l). e 2 {V h u(z ia )) 2 . (13f) 

The curve T will grow locally at x(0), if the curve moves 
in direction —r(0). In this case Tt a n(0) = £t(0). r(0) < 0. 
Therefore, 

a< |f(0).e 1 |(V 2 u(i O ’ 1 )) 2 +|f(0).e 2 |(V^’ 2 ))| 2 (14) 

has to be satisfied such that the curve length increases. 
For the exemplary case r(0) = ei, the condition reduces 
to 0, i.e. to the condition from the introductory exam¬ 
ple. 

The curve T(t) will grow at x(l), if the curve moves 
in direction r( 1) leading to Vtan(l) > 0. Therefore, the 
inequality 

a < 1^(1) • 61 |(V^u(i ' 1 ’ 1 )) 2 + |r(l). e 2 |(V^u(7 1 ’ 2 )) 2 (15) 


has to be satisfied. 

Since the term cr|T| in the energy 0 penalizes the 
length of the curve, a curve can only grow in direction 
—r(0) or r( 1), if the derivative terms (V^tx) 2 , i = 1,2, 
are large compared to a. 

The scheme (13) describes the motion of the curve. For 
Nc curves Ti,..., Tjsf c , we can solve (13) for each curve. 
For a closed curve only the normal velocity (13a) with 
the relation 


([13b]) needs to be considered, on noting that 
x^(0) = Xi(l). In the case of triple junctions and intersec¬ 
tions with the image boundary, additional conditions for 
the involved endpoints have to be considered. If triple 
junctions occur, the evolution equations for the corre¬ 
sponding three curves which meet at the junction are 


For computing the position of the evolving curves T nu¬ 
merically, we consider a decomposition of the interval 
[0,1] of the form 0 = q l 0 < q\ < ... < q l N . = 1 , for 
i = l,..., Nc • If is a closed curve, we make use of the 
periodicity Ni = 0, Ni + 1 = 1 , — 1 = Ni — 1 , etc. 

Further, let 0 = to < t\ < ... < = T be a 

partitioning of the time interval [0, T] with time steps 
A t m := t m+ 1 — t m , m = 0,. .., M — 1 . Smooth curves 
r i (t m ), 7 = 1 ,...,Ac, m = 0 ,...,M are replaced by 
polygonal curves T™ given by nodes which are ap¬ 
proximations of x^q^tm)- Further, let be an approx¬ 
imation of ^i(qj^t m ). The derivative terms with respect 
to time are replaced by difference quotients of the form 

* -N- (x™ +1 - X%) . (16) 

Let h ™._t = XJj - XJj-i, i = l,...,N c ,j = l,...,N i , 

be the distance between two neighboring nodes. For each 
curve, we define a discete normal vector field VJ 12 by 


Vi 


f vm vm 
\ X i,3 l) 






see also J6j7,9 . Here, _L denotes the anti-clockwise rota¬ 
tion of a vector by tt/2. Further, we define the following 
weighted approximating normal vector at X™j by 


h m 


* %" v •— 
U i, j 


l ,3 2 Z ’-7 2 


+ h™ t iZ n .. 1 

l ,3~\~ 2 l i3 ~2 


( vm vm \ J 




i,3 +1 


M- 




for 7 = 1 ,..., Ni if dT™ = 0 and for j = 1 ,..., Ni — 1 if 
dT™ 7 ^ 0. In the latter case, we set 


m \ JL 


q,o 


(X™-X™ 0 ) 
h n 


“i,Ni := V i 


fvm _ vm \ - 
n = 1) 

, 2 Tfl 

n i,Ni-± 
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The external term F is approximated by 

\u 0 {X^)-u{X^+a^)) 2 
-(uoiX^-uiX^-a ^)) 2 


Ki : = A 


with a small real number a > 0 if j is not the index of a 
free endpoint, and we set F-j = 0 else. 


The equation for the normal velocity (13a) is approxi¬ 
mated by 




rrK rn + 1 
3 aK iJ 


+ F\ 


i,3' 


(17a) 


Thus, for computing T^ +1 , we use the previous curve T™ 
for the external term F™ and for the weighted normal 

wT 






For an approximation of (13b), we need to define an 
approximation of (xi) ss (qp £ m+ i). For that, we make use 
of difference quotients of the form 

o 

a h,m -y-m+l __ _ 

2 i,j ■ h m 

l ,3 2 2 

-(X^-X^/h™^), 

for i = 1 ,..., N c and j = 1 ,...,W, if dT? = 0, 
and j = — 1, else. In case of equal spatial 

step sizes h m . x = h m .. 1 =: FF , the term reduces to 

- 2 + X™- +1 ) / ((h™) 2 ), see also |5], where we 
also defined and used these difference quotients. 


The equation (13b) is now approximated by 


m+l.-jm _ a h,m C?m+1 

H,j u i.i ~ A 2 A i.i > 




(17b) 


for i = 1 ,..., TVc, j = 1 ,..., N{ in case of closed curves 
and j “ 1 ,..., Ni — 1 in case of open curves. 

In case of open curves, additional equations for the 
endpoints are needed. The case of triple junctions and 
boundary intersection points is described in [9 . For 
curves with free endpoints we introduce the tangential 
vectors = (X% - X™ 0 )/\ h , and f™ Ni = (X™ Ni - 


Xr Ni _i)/h iN ._i. The equations (13c) and (13d) are ap¬ 
proximated by 

=a-\T™ 0 .e 1 \(X 2 h u(^)) 2 - |ry 0 .e 2 |(V^’ 2 )) 2 , 

(17c) 


and 

^ ( -ym+1 

A t m v W 


x 


m \ =?rr 

Xi) ■ T i , 


Ni 


= -v + K Ni - eiKV^f 1 ’ 1 )) 2 + \r™ Ni . e 2 |(V^(#’2))2. 

(17d) 


Similarly, discrete versions of (13e) and (13f) can be 


The scheme is a numerical approximation of the 
scheme (13), where the parametric curves are replaced by 


polygonal curves, and the smooth functions Xi and are 
replaced by continuous functions uniquely given by their 
values at the nodes q 1 ^ i = 1 , ..., Nc , j = 0 ,..., N{. 

The discrete scheme can be rewritten to a linear system 
with a sparse system matrix, similar as presented in [9 , 
and can be solved with a fast direct solver like for example 


the UMFPACK algorithm 19 


3.2 Numerical Solution of the Denoising 
Problem 

For computing a numerical solution u h for the piece- 
wise smooth, denoised version u of uo, we consider for 
X x , N y e N the discrete set 

n h := {( ih,jh ) : i = 0 ,..., N x , j = 0 ...,N y }, 

where N x and N y are the number of pixels in x- and 
^/-direction. We define for i = 1,..., N x , j = 1,..., N y 

h 2 , if [(i-i)h,ih]x{j}nr% = 9, 
A x {i,j)-{ Vi 0 e { 1 ,..., Nc}, 

0 , else, 


h 2 , if (0 x [(j-l)h,jh]nT% = t 
A y( i d) = { Vi 0 e { 1 ,..., Nc}, 

0 , else. 


Fixing the set of curves T, we consider the following 
discrete energy: 


N X Ny 

-E'discr (N ) ^ [ A x (i,j) 

i =1 3 =1 


u ij ~ u i-i,j 


+ Ay(i,j) 


,j — 1 

h 


N X Ny 

+ a££ h 2 (uo(ih,jh) — u^) 2 , (18) 

*=o j=o 

which is a discrete analogue of 
f n \ p (|| V ^|| 2 dx + f Q A(u 0 - u ) 2 ) dx. Here u\ - ap¬ 
proximates u at the node ( ih,jh ). The piecewise 
continuous function u h is uniquely given by its value at 
the points in Q h . 

By setting the terms A x (i,j) or A y (i,j ) to zero at 
points where the line [(i — 1 )h,jh), (ih,jh)] or [(i/i, (j — 
1 )h),(ih,jh)] intersects with one of the curves, we ap¬ 
proximate the integral over the set Q\T. 


Taking the derivative of the right hand side of (18) with 


stated using an d ^Tn- as discrete normal vectors. 


respect to u^ - and setting the resulting term to zero, leads 
to a linear system. The corresponding system matrix is 
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sparse since each node ( ih,jh ) G fl h is only coupled to 
a few neighboring nodes. The resulting linear system 
can be solved with a fast direct or iterative solver by 
employing the sparse matrix structure. 

Considering h 0, we obtain in the limit Xu . V = 0 at 
the curves belonging to T, and Xu . tiqq = 0 at the image 
boundary dQ, where tiqq is a normal vector field at dCl. 
For details, we refer to [9]. Consequently, we obtain an 
edge preserving image smoothing if T matches with the 
edges in the given image. 


2. Compute the external terms F-j by using the solu¬ 
tion u h of step [l] Compute X m+1 by solving the 
linear equation derived from the scheme 

3. Check whether topology changes occur. If so, exe¬ 
cute the topology change. 

A segmentation of the image is given by the final set of 
curves V M . An image restoration is given by the image 
approximation u h from the time step 


3.3 Topology Changes 

During the evolution of curves, topology changes can oc¬ 
cur, since the edge set in the image and the boundaries 
of objects are not known in advance. Therefore, curves 
can split into two or more subcurves, curves can merge 
to one single curve, triple junctions and new curves may 
occur and curves can intersect with the image bound¬ 
ary such that new boundary nodes emerge. Further, a 
curve needs to be deleted if its length becomes too small. 
In 9 , we extended the idea of [5,26 , and described a 
method to detect topology changes of curves efficiently. 
The main idea is the use of an artificial background grid 
which covers the entire image domain D. We consider 
successively all nodes and mark a grid element with 
(i,j) if is the first node located in this array. If a 
grid element is already marked with and the nodes 

X™j and X™ ^ are not neighbor nodes, a topology change 
likely occurs close to this pair. Details on this method 
for curves without free endpoints are given in [9 . 

In principle, topology changes involving curves with 
free endpoints can be detected similarly by using such a 
background grid. In addition to the topology changes 
listed above (splitting, merging, emergence of triple 
junctions and boundary intersection points), topology 
changes involving the free endpoints can occur: If two free 
endpoints of one curve are located in one square of the 
background grid, an open contour becomes a closed con¬ 
tour. If two free endpoints of two different curves meet, 
the two curves merge to one single curve, and the former 
free endpoints become inner nodes of the new curve. If a 
free endpoint and an inner point of a curve meet, a triple 
junction is created. 


3.5 Modifications 

Step [2] of the algorithm above can additionally be split 
in two sub-steps: First, we fix the free endpoints and 
we let the inner nodes of the curve evolve. Then, we let 
the endpoints evolve according to the above presented 
discrete scheme. 

The main effort of this method compared to the Chan- 
Vese method for interface curves is that we have to solve 
a two-dimensional diffusion equation (bulk equation) sev¬ 
eral times during the segmentation. In the experiments 
described in the next section, we perform 10 steps of 
curve evolution followed by a solution of the bulk equa¬ 
tion. Having computed u h , we use it for the next 10 curve 
evolution steps. 

As an alternative, we can start the segmentation using 
interface curves and the image segmentation method de¬ 
scribed in 19] (based on the Chan-Vese method [l5]) with 
piecewise constant approximations. As a postprocessing 
step, we can consider the derivatives of the image func¬ 
tion in normal direction at the final curves (or the jump 
of the image function across the curves). We replace in¬ 
terface curves by curves with free endpoints if the deriva¬ 
tives in normal direction are locally very small. For that, 
we delete those parts of a curve where the derivative is 
small which results in curves with free endpoints. Next, 
we compute some steps of the segmentation method with 
free endpoints to obtain the final contours. 

Topology changes occur only in rare cases when using 
a postprocessing evolution of curves with free endpoints. 
In most situations, topology changes are already detected 
in the previous evolution. 

4 Results and Discussion 


3.4 Summary of the Algorithm 

We propose the following algorithm for image segmenta¬ 
tion and image restoration with parametric contours with 
possible free endpoints: 

Given a set of polygonal curves T° = (T®,... ,T% c ) 
and X° = (X°,... ,X% c ) with X?([0,1]) = T?, perform 
the following steps for m = 0,1,..., M — 1: 

1. Compute a denoised image approximation u h by 
minimizing ( [l8| ) (solve the corresponding sparse lin¬ 
ear system). 


The method for image segmentation and restoration pre¬ 
sented in sections [2] and [3] are applied on some exem¬ 
plary test images. For all experiments presented in this 
section, we use constant time steps sizes A t m = At, 
m = 0,..., M — 1. 

In the first experiment, we consider an example where 
a contour with two free endpoints evolves in the image 
domain and detects an edge. Figure [3] presents the results 
of image segmentation and denoising. It can be observed 
that the image is not smoothed out across the curve T. 
Further a growth of the curve in tangential direction can 
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Figure 3: Example image showing a contour with two 
free endpoints. Original image and evolving contours ( 1 st 
row) and denoised image ( 2 nd row) for m = 1 , 1000 , 6000 
using At = 0.032, a = 2 e — 5 and A = 0 . 002 . 


Figure 4: Example image showing a contour with one 
free endpoint and one boundary intersection point, see 
also 32 . Evolving contours for m = 1,500,3000 using 
At = 0 . 001 , a = 2 e — 5 and A = 0 . 002 . 



be observed. The growth stops when the inequalities (14) 


and (15) become equalities. This depends on the absolute 


values of the difference quotients \S7 l h u\, i = 1 , 2 , and 
the weighting parameter a. The image approximation 
u attends values in [0,1]. In this image, differences of 
the form u(x + hei) — u(x) are typically of magnitude 
10 -2 . Since Q = [1,300] x [1,300] and h = 1 , \S7\u\ 2 is 
of magnitude 10 -4 . Therefore, we have to choose a small 
value for the weighting parameter cr, here, we choose a = 
2e —5. If we used a normalized image domain D = [ 0 , 1 ] x 
[0,1], the pixel grid would have a grid size of h = 1/300 
and h 2 = 1/90000. In this case, we would choose a weight 
a of magnitude 1 . 

In a second experiment, we study a crack tip problem 


which has also been considered in 32 . The image func¬ 
tion is given by 


uq(x ) = a\J r(x) sin(#(x)/ 2 ) + 5, 


(19) 


where r{x) > 0, Q(x) G (—7r, 7r] are polar coordinates with 
r = 0 corresponding to the image center, and a, b G M are 
constants such that uq attends values in [0,1]. 

Figure [4] shows the evolution of a contour with one 
free endpoint. The second endpoint belongs to the image 
boundary. At time step m = 3000, the free endpoint is 
located at the image center and the curve matches with 
the edge in the image. As discussed above (see also Equa¬ 
tions © and ©), the parameter cr, which weights the 
length term, needs to be chosen small enough such that 
the curve can extend. If a is fixed, the absolute value of 
the difference quotients must be large enough such that 
the length of the curve increases. In this example, the 
edge is a horizontal line and the position where the curve 
stops depends on the value of V 2 rq i.e. on the difference 
quotient in ^/-direction. 

We rerun the example using a = 0.002 and a = 0.01 
instead of a = 2e — 5. Figure [5] shows the results at time 


Figure 5: Dependency on the weighting parameter using 
a = 0.002 (left) and <r = 0.01 (right), At = 0.001 and 
A = 0.002 for m = 3000. If a too large weight is chosen 
for the length term in (|7|), the curve does not reach the 
image center. 


step m = 3000. In both cases, the free endpoint does 
not reach the center of the image since the value of a has 
been set larger. The growth of the curve already stops 
at larger values of V|rq recall conditions (14) and (15)). 
We even let the curve evolve until time step m = 5000, 
but we observed no significant motion between m = 3000 
and m = 5000. 


In another experiment, which demonstrates the evo¬ 
lution of curves with free endpoints, we first apply the 
parametric method of 9] to the Chan-Vese problem 
using interface-curves. This means, that we first seg¬ 
ment a given image in regions separated by interface 
curves, see Figure [6j We start with one large initial curve 
which splits up in two sub-curves. This example thus also 
demonstrates the handling of a topology change. 

In a postprocessing step, we delete those nodes where 
the jump of uo across the curve is smaller than a given 
tolerance of tol = 0.1. This results in one closed curve, 
where no points are deleted (blue curve in Figure [7]), and 
in one curve with two free endpoints (red curve). Fig¬ 
ure [7] shows the results of a postprocessing evolution of 
the curve. This example shows that our methods for im¬ 
age segmentation and denoising can be applied also on 
images with both open and closed edges. 

Table [l] shows the values of the discrete Mumford-Shah 
energy (|7|) for the last step of the Chan-Vese piecewise 
constant segmentation with closed region boundaries (cf. 
Figure [ 6 j m = 1200) and for the initial and final step 
of the postprocessing with one open boundary (cp. Fig- 
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Figure 8: Postprocessing evolution with a contour with 
two free endpoints using tol = 0.5 to obtain the initial 
contour. Original image and contours (left, center) for 
m = 1, 600 using At = 0.05, a = 2e — 5, A = 2e — 4 and 
denoised image for m = 600 (right). 



Figure 6: Image segmentation result using Chan-Vese, 
interface curves and a piecewise constant image approx¬ 
imation. Original image and contours (first row) and 
piecewise constant approximation (second row) for m = 
1,1050,1200. 
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Figure 7: Postprocessing evolution with a contour with 
two free endpoints using tol = 0.1 to obtain the initial 
contour. Original image and contours (left, center) for 
rn = 1,400 using At = 0.05, a = 2e — 5, A = 2e — 4 and 
denoised image for m = 400 (right). 


Table 1: Comparison of discrete Mumford-Shah Energy 


Processing 

Method 

Step Nr 

Discrete 
Mumford- 
Shah En¬ 
ergy 

Chan-Vese, 

1200 (fi¬ 

22364.94 

piecewise con¬ 
stant 

nal) 


Postprocessing, 
free endpoints 

1 (start) 

23162.04 

Postprocessing, 

400 

18284.36 

free endpoints 

(final) 



ure [7J m = 1 and m = 400). Note, that the absolute 
values are large, since the image consists of 90000 pixels 
and the size of each pixel is 1 x 1 . The average contribu¬ 
tion of each pixel to the energy is < 1; however, it sums 
up to a value of magnitude 2 • 10 4 . 

From Table [l] we can observe that the energy even 


slightly increases from the last step of the Chan-Vese 
piecewise constant method to the first step of the post¬ 
processing evolution. Deleting part of the curve does not 
decrease the energy in this example. However, at the end 
of the postprocessing evolution with one open contour, we 
obtain a decrease of the energy. The energy at m = 400 
of the postprocessing is 81.75% of the energy of the last 
step of the piecewise constant segmentation with closed 
boundaries. Therefore, if we delete part of the curve and 
if we let the curve with free endpoints evolve again, we 
will obtain a final curve such that the corresponding dis¬ 
crete Mumford-Shah energy 0 is reduced compared to 
the Chan-Vese piecewise constant result. 

Next, we investigate the influence of the tolerance value 
to/, which is used for the deletion of some nodes of the 
curve. Note, that the image function attends values in 
[0,1] where 0 corresponds to black and 1 corresponds to 
white color. The exact value of the tolerance tol influ¬ 
ences only the start curve of the second curve evolution. 
We repeat the postprocessing evolution and use tol = 0.5 
as tolerance resulting in a different initial curve. 

Figure [8] shows the postprocessing evolution of the 
curve. Of course, since the initial contour in Figure [8] 
(left) is smaller compared to the initial contour in Fig¬ 
ure [7| more iteration steps are needed to obtain the final 
contour. 

The final result is independent on the exact initial 
curve as long as it is of the same type, i.e. open with 
free endpoints; not closed or not fully deleted. In our ex¬ 
ample the largest jump of uq across the final red curve of 
the first evolution (cf. Figure [6j right sub-figures) is 0.61, 
the smallest jump is 0.05. The large difference between 
the largest and smallest jump can be used as an indicator 
to replace the interface-curve by a curve with two free 
endpoints. (On the contrary, the jump across the blue 
curve is constant in this example.) As tolerance value tol 
any value larger than 0.05 and smaller than 0.61 could be 
chosen. For tol < 0.05 no node point would be deleted 
resulting in an unchanged curve. For tol > 0.61 all nodes 
and therefore the entire curve would be deleted. All val¬ 
ues between the two thresholds can be theoretically used. 
Therefore, in this example, the final result is independent 
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Figure 9: Medical image segmentation using Chan-Vese, 
interface curves and a piecewise constant image ap¬ 
proximation. Original image and contours for m = 
1,400,1000,2500. Image courtesy: Dr. Declan O’Regan 
and the Robert Steiner MR Unit, MRC Clinical Sciences 
Centre, Imperial College London. 



Figure 10: Postprocessing segmentation with free end¬ 
points. Original image and contours (far left, left, right) 
for m = 1,100,500 using At = 0.008, a = 3.3e — 5, 
A = 0.0167 and denoised image (far right) for m = 500. 
Image courtesy: Dr. Declan O’Regan and the Robert 
Steiner MR Unit, MRC Clinical Sciences Centre, Impe¬ 
rial College London. 


on the exact value of tol as long it is in (0.05, 0.61). 

Next, we demonstrate an example where a real, med¬ 
ical image is processed. Figure [9] and Figure [lO] show 
an excerpt of a medical image and the result of an edge 
detection. We first use the Chan-Vese algorithm with 
piecewise constant image approximation for segmenting 
the image, see Figure [9j In this first segmentation step, 
also topology changes occur. The initial closed curve 
touches twice the image boundary and splits up in two 
open curves each with two boundary intersection points. 
After the preceding segmentation, a part of the red curve 
is deleted (using a tolerance of 0.1 for the jump across the 
curve) resulting in a curve with free endpoints. Figure [l()| 
shows the result of the postprocessing evolution. Small 
tangential motions of the free endpoints can be observed. 

Finally, we study an example where several topology 
changes occur. Figure |TT] shows an example where we 
start with many small initial curves. A similar image is 
also considered in 32 . Many of the small initial lines 


shrink and are deleted when their curve length becomes 
too small. Additonal topology changes occur: Near the 
upper left corner of the image, two curves merge at their 
free endpoints to one curve. Further, three free endpoints 
become boundary intersection points, and a triple junc¬ 
tion emerges when a free endpoint meets another curve 



Figure 11: Image segmentation and contour detection 
with topology changes. The final segmentation contains 
two free endpoints, one triple junctions and three bound¬ 
ary intersection points. Original image and contours for 
m — 1,100, 250,400, 750,1200 using At = 1, a = le — 4, 
A = le — 3. 


at an inner node. The topology changes are detected as 
described in Section 3.3 and 9 . 


An advantage of the parametric method is that we 
can easily handle non-interface curves and complex curve 
networks including triple junctions. The curve evolution 
scheme is very similar to the scheme presented in [9 for 
interface-curves. Instead of computing the mean value of 
the image function in regions, we have to solve a diffusion 
bulk equation. Additional to the motion of the curve in 
normal direction, free endpoints can move in tangential 
direction. 


There are alternatives to parametric methods to de- 

is very 


scribe an evolving curve. The level set method 30 


popular for image processing applications and in partic¬ 


ular for active contours methods, see e.g. 25 , 11 , 23 


15 , 37 , 35 to mention a few. In level set methods, a 


hypersurface is embedded as the zero level set of a func¬ 
tion defined on the image domain D. With level set tech¬ 
niques, free endpoints however cannot be handled with 
one single level set function: Since level set methods em¬ 
bed a curve as zero level set of a function 4> : U M, the 
curve is an interface between two regions {4> > 0} and 
< 0}. Therefore level sets are always closed or meet 
the image boundary at their endpoints. Non-interface 
curves can be handled by using two level set functions <f> 
and T and by using artificial regions, see [36 . A curve 
with free endpoints can then be represented by the inter¬ 
face between the artificial regions {4> > 0} fi {T > 0} and 
> 0} D {T < 0}, for example. 


Using our direct, parametric approach it is not neces¬ 
sary to introduce artificial regions. Further our method 
is very efficient, since the curve evolution is only a one¬ 
dimensional problem. 
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5 Conclusion 

We proposed a new parametric approach for active con¬ 
tours with free endpoints. The image segmentation and 
denoising method presented in this article is based on 
a discrete version of the Mumford and Shah functional. 
For curves with free endpoints a flow in normal direc¬ 
tion and a flow of the endpoints in tangential direction 
attracts the curves to the edges in the image. With the 
presented approach, we can handle both open and closed 
curves. The method is also suitable to be employed as 
postprocessing step to improve the result of a previous 
Chan-Vese like segmentation with interface curves and 
piecewise constant image approximations. 
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